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Abstract 


This report documents the code AMR1D, which is currently posted 
on the World Wide Web [ http://sdcd.gsfc.nasa.gov/ESS/exchange/co 
ntrib/de-fainchtein/adaptive_mesh_refinement.html ]. AMR1D is a 
one dimensional finite element fluid-dynamics solver, capable of adap- 
tive mesh refinement (AMR). It was written as an instructional tool 
for AMR on unstructured mesh codes. It is meant to illustrate the 
minimum requirements for AMR on more than one dimension. For 
that purpose, it uses the same type of data structure that would be 
necessary on a two dimensional AMR code (loosely following the al- 
gorithm described by Lohner [1]). 




1 INTRODUCTION 

The basic idea of AMR is to provide high resolution to a simulation, only 
WHEN and WHERE it is needed. For that purpose, a mechanism to eval- 
uate the need for refinement “on the fly” is necessary. AMR1D uses the 
“Error Indicator” function of Lohner [1], This “Error Indicator” is an upper 
bound estimate of the L 2 norm [2] of the discretization error, which has been 
normalized and filtered for noise. 

Unstructured codes are especially well suited for adaptive mesh refine- 
ment. They allow for the total disengagement of the mesh refinement pro- 
cess from the solver algorithm. When AMR executes it delivers a (modified) 
unstructured mesh, with all the appropriate mesh variables, as well as values 
for all the unknowns at each mesh point. Because the fluid solver is already 
equipped to handle unstructured meshes, this mesh information is sufficient 
for the fluid solver to advance the solution. The fluid solver requires no 
additional information pertaining to the refinement history of the mesh. 

AMR1D is designed to return the original mesh when and where high 
resolution is no longer needed, through adaptive DE-REFINEMENT. Thus, 
when a portion of the mesh is de-refined, AMR1D simply undoes a previous 
mesh refinement. A “history array” that keeps track of the refinement history 
of the mesh is maintained for that purpose. 

AMR1D allows for multiple mesh refinements. The maximum number of 
refinements that a given element may undergo is determined by a user-defined 
parameter. However, a given element may only be refined or de-refined once 
each time AMR is invoked. 

Once the elements requiring refinement are selected, a layer of elements 
around them are also selected for refinement. This layer acts as a buffer to 
maintain accuracy at high resolution regions. 

The contents of this report are organized as follows: Section 2 provides 
the general algorithm of AMR1D. Section 3 explains briefly the flow-solver 
algorithm used. Section 4 explains in some detail the AMR algorithm, and 
the subroutines that compose it. Section 5 is a description of the variables 
and parameters used in the code. Section 6 contains the Fortran code for 
AMR1D. 

Figure 1 shows results produced by running AMR1D with initial condi- 
tions set for the Sod Shock Tube Problem [3]. 
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Sod Shock Tube Problem (DENSITY) 




Figure 1: Results produced by AMR1D for the Sod Shock Tube problem. The 
first plot shows the density. Below it, the refinement level of each element is 
plotted. 
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2 THE GENERAL ALGORITHM 


AMR1D consists of a loop over the time steps, preceded by some initialization 
routines. 

INITIALIZATION 

• call INIT 

Define all simulation and AMR parameters, as well as initial conditions. 

• call GEOM 

Compute the mesh geometry and connectivity array values (e.g. ele- 
ment lengths). 

• call AMR ROUTINES 

Refine a few times before beginning the simulation. 

MAIN LOOP OVER TIME STEPS 

• If NSTEP = INTEGER x NREF => call AMR ROUTINES 

• Advance the solution with the FLOW-SOLVER ROUTINES 

• Apply Boundary Conditions 

• If NSTEP = INTEGER x NWRITE =► write output 

3 THE FLOW-SOLVER ALGORITHM 

This is a conservative algorithm. It solves a set of three equations of the 
form 

du/dt + df/dx = 0 (1) 

by advancing in time the three conserved variables: mass density (u=rho), 
linear momentum (u=rhov), and specific total energy (u=rhoE). The appro- 
priate fluxes (f) at the points are computed in subroutine FLUX1D. They 
are the fluxes corresponding to the Euler equations. 
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ie 


•o — # # — a 

ipl ip2 

intmat ( ie , 1 ) =ipl 
intmat ( ie, 2 ) =ip2 

Figure 2: Schematics of the grid connectivity. The two nodes of element ie 
are the points numbered ipl and ip2, respectively. 


The time step is also computed in FLUX1D using a Courant condition: 
dt = COUR * min[RLEN/v m } (2) 

where, 

Vm = V s + M, (3) 

v s is the local speed of sound and v is the local velocity. The minimum is 
computed over each element, COUR is a user-defined parameter (between 0 
and 1), and RLEN is the length of the element. 

The increments (du) for each time step are computed in subroutine 
DELTAU, using a Galerkin finite element scheme. The elements are La- 
grangian [4] and have two nodes (points) and piecewise linear interpolation 
functions. Figure 2 illustrates the grid connectivity between elements and 
points. The scheme is first order in time and space. First order numeri- 
cal diffusion is added for stability, in a manner similar to that used in the 
Rusanov scheme [5, 6]. 

The flow-solver algorithm is modeled after UNST, a simple two dimen- 
sional unstructured fluid solver written by S.Zalesak 1 

4 THE AMR ALGORITHM 

The refinement algorithm consists of three steps, performed by calling three 
subroutines: 

1 Private Communication 
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• subroutine GRERROR 


• subroutine UNREF 

• subroutine REFINE 

Each element can only be refined or de-refined once, during a given pass 
through the refinement routines. After refinement, subroutine GEOM is 
called to re-compute all the geometry arrays for the “new” mesh. 

SUBROUTINE GRERROR 

Subroutine GRERROR determines which elements should be refined or de- 
refined. The Error Indicator function (ERRORE) is computed and assigned 
to each element (ie). An element with a value of ERRORE larger than 
the user defined threshold for refinement (CTORE), is marked for refine- 
ment: LREFE(ie) = l. Similarly, an element with a value of ERRORE smaller 
than a user defined threshold for de- refinement (CTODE), is marked for de- 
refinement: LDERE(ie)=l. 

Next, a user-defined number of buffer layers of refined elements (NBUFF) 
are added around each element originally marked for refinement. 

Finally, a check is performed to avoid “jumps” in the level of refinement 
on neighboring elements. This step is only necessary on one-dimensional 
grids. Connectivity rules on two and three dimensional grids automatically 
avoid discontinuity on the level of refinement of neighboring elements. 

SUBROUTINE UNREF 

Subroutine UNREF de-refines elements marked for de-refinement in GR- 
ERROR. Mesh de-refinement follows two rules: 

• Elements are de-refined, as opposed to coarsened. That is, a pair of 
sibling elements marked for de-refinement are not arbitrarily coarsened, 
but are instead replaced by the original “parent” element. (Note that 
the initial mesh is the coarsest mesh allowed). This rule makes it 
necessary to maintain a refinement history element array: MRHIST. 

• In order for de-refinement to take place, both siblings must be marked 
for de-refinement. 
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The elements marked for de-refinement in subroutine GRERROR are checked 
for consistency with these rules before being de-refined. The actual de- 
refinement proceeds as follows, 

1. The de-refined elements are de-activated: MRHIST(ie_child,4)=0 

2. The parent element is activated: MRHIST(ie_parent,4)=l 

3. The “extra” point (jp) previously created to produce the sibling el- 
ements out of the parent element (c.f. SUBROUTINE REFINE) is 
de-activated: IPACT(jp)=0 

SUBROUTINE REFINE 

Element refinement consists of dividing the parent element into two children 
elements of equal length, sharing a newly created point. 

Subroutine REFINE refines some or all of the elements marked for re- 
finement in subroutine GRERROR. AMR1D allows for only a (user-defined) 
given maximum number of refinements of each original element(NMXRE). 
Thus, subroutine REFINE checks every element marked for refinement in 
GRERROR, and refines only those whose level of refinement will not exceed 
NMXRE. 

The refinement process is fairly straightforward: 

1. Add (activate) a new point at the end of the point list. (The point will 
be at the center of the element ie marked for refinement). 

2. By linear interpolation, compute the new coordinate and unknowns: 
density, momentum and specific total energy. 

3. Modify the connectivity array (INTMAT) to add the two new elements. 

4. Disable the parent element and enable the children elements. 

5. Update the refinement history element array (MRHIST). 

6. Update NPOIN and NELEM. 
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5 DEFINITIONS 


FLOW SOLVER VARIABLES 


NPOIN 

NELEM 

INTMAT 


XP(ip) 

RHO(ip) 

RHOV(ip) 

RHOE(ip) 

CE(ie,in) 

RLEN(ie) 

RMATM(ip) 


FLUX(ip,l:3) 

P(ip) 

V(ip) 

PERIOD 


= the number of points 
= the number of elements 
= the “element list” 

( INTMAT(ie,l) and INTMAT(ie,2) are the first and second point 
numbers that comprise element ie) 

= the coordinate of point ip 
= the density at point ip 
= the momentum at point ip 

= the specific total (internal plus kinetic) energy at point ip 
= coefficients for evaluating x first derivatives on an element 
= the length of element ie 

= the ip diagonal element of the inverse lumped mass matrix 
( = 2./ [sum of lengths of elements sharing point ip] ) 

= flux functions f for point ip 
= pressure at point ip. 

= velocity at point ip 

= mesh period length (when PERIOD=0, boundaries are “hard walls”). 


ADAPTIVE MESH REFINEMENT VARIABLES 


ERRORE(ie) 

= “Error indicator” 


MRHIST(ie,l) 

= parent element number. 

(default=0) 

2) 

= first child element number. 

(default=0) 

3) 

= second child element number. 

(default=0) 

4) 

= level of refinement. 

(default=0) 

5) 

= element active? 

( 1 = yes ) 

LREFE(ie) =1 

=>■ refine 


=0 

=*> do not refine 


LDERE(ie) =1 

=> de-refine 


=0 

=$■ do not de-refine 


IPACT(ip) =1 

=> point ip is active 


=0 

=> point ip is not active (deleted) 
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CTORE = threshold for refinement 
CTODE = threshold for de-reninement 
EPSIL = noise tolerance parameter 
NBUFF = No. of buffer layers 

6 THE FORTRAN CODE 

program AMR ID 
c 

c written by: Rosalinda de Fainchtein 

c 

c 

parameter(mpoin=40000,melem=80000,period=0 .0) 
integer intmat(melem,2) 

real xp(mpoin) ,ce(melem,2) .rmatm(mpoin) ,rlen(melem) 
real rho(mpoin) ,rhov(mpoin) .rhoE(mpoin) 
real p(mpoin) ,v(mpoin) ,f lux(mpoin,3) 
real diff (melem,3) 
real vmax(mpoin) 
real du(melem,3), dup(mpoin,3) 
c 

integer mrhist(melem,5) , lrefe(melem) , ldere(melem) 
integer ipact(mpoin) 
integer iptemp(mpoin) ,ietemp(melem) 
integer l(mpoin) 

real tempea(melem) ,tempec(melem) .errore(melem) 

real temppa(mpoin) .temppb(mpoin) ,temppc(mpoin) .errorp(mpoin) 


c 

c Get initial conditions and parameters 

c 

call init (gamma , cour , ntime ,di , ctore , ctode , epsil , nrmax , 
& nbuff ,ntref ,npo in, nelem,mpoin,melem, period, 

& intmat, mrhist, lrefe, ldere, ipact, 

& xp, rho, p, v, rhov, rhoE, 

& ibdry_l, ibdry_r ) 
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gamma l=gamma-l . 
c 

c Calculate element lengths, inverse mass matrix, and 

c coefficients for evaluating x derivatives on the elements, 

c 

call geom(mpoin, melem, npoin, nelem, intmat, xp, ce, rmatm, 

& rlen, period, mrhist, ipact ) 

c 

c Refine a few times before starting 

c 

do iloop=l,5 
c 

c Refinement Routines 

c 

c 

c 

c Compute the error function at each element and mark elements 

c that require refinement or de-refinement 

c 

call grerror (nelem, npoin, melem, mpoin, rho, 

& tempea, tempec, temppa, temppb, temppc , 

& epsil, ctore, ctode, nbuff ,errorp, 

& errore, mrhist, ipact, 

& lrefe, ldere, ce, rlen, intmat, 

& ietemp,iptemp,ibdry_l,ibdry_r ) 

c 

c Coarsen elements marked for derefinement 

c 

call unref(melem, mpoin, nelem, npoin, intmat, 

& mrhist, ldere, ipact ) 

c 

c Refine elements marked for refinement 

c 

call ref ine (melem, mpoin, nelem, npoin, intmat, nrmax, 

& mrhist, xp, rho, rhov, rhoE, lrefe, ipact) 

c 
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c Re-compute element lengths, inverse mass matrix and coefficients 

c for computing x derivatives within each element 

c 

call geom(mpoin, melem, npoin, nelem, intmat , xp, ce, rmatm, 

& rlen, period, mrhist, ipact ) 

c 

c End of Refinement Routines 

c 

c 

end do 
c 

c BEGIN THE MAIN LOOP OVER TIME STEPS 

c =================================== 

c 

do 1000 i=l,ntime 
c 

c Refinement Routines 

c 

c 

if ( mod(i,ntref) .eq. 0 ) then 
c 

c Compute the error function at each element and mark elements 

c that require refinement or de-refinement 

c 

call grerror (nelem, npoin, melem, mpoin, rho, 

& tempea, tempec, temppa, temppb, temppc , 

& epsil, ctore, ctode, nbuff, errorp, 

& errore, mrhist, ipact, 

& lrefe, ldere, ce, rlen, intmat, 

& ietemp, iptemp ) 

c 

c Coarsen elements marked for derefinement 

c 

call unref (melem, mpoin, nelem, npoin, intmat, 

& mrhist, ldere, ipact ) 

c 

c Refine elements marked for refinement 
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c 

call ref ine (melem, mpoin, nelem, npoin, intmat, nrmax, 

& mrhist, xp, rho, rhov, rhoE, lrefe, ipact) 

c 

c Re-compute element lengths, inverse mass matrix and coefficients 

c for computing x derivatives within each element 

c 

call geom(mpoin, melem, npoin, nelem, intmat, xp, ce, rmatm, 

& rlen, period, mrhist, ipact ) 

end if 
c 

c End of Refinement Routines 


c 

c 

c 


c Compute fluxes at the points and allowed time step. 

c 

call fluxld (mpoin, melem, npoin, nelem, gamma, cour, 

1 intmat, rho, rhov, rhoE, p, v, 

2 rlen, vmax, flux, dt, mrhist, ipact ) 
c 

c Compute changes in conserved variables on elements, and 

c scatter-add them to the point array dup 

c 

call deltau (mpoin, melem, npoin, nelem, intmat, du, dup, flux, dt , 

& ce.rlen, cour, di, rho, rhov, rhoE, dif f , vmax, mrhist ) 

c 

c Hard-Wall boundary conditions 

c 

if (period .eq. 0. ) then 
dup(ibdry_l ,2)=0 . 
dup ( ibdry.r , 2) =0 . 


end if 
c 

c Advance the solution. 

c 

do ip=l, npoin 
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if ( ipact(ip) .eq. 1 ) then 
rho(ip) = rho(ip) - dup(ip,l)*rmatm(ip)*dt 
rhov(ip) = rhov(ip) - dup(ip,2)*rmatm(ip)*dt 
rhoE(ip) = rhoE(ip) - dup(ip,3)*rmatm(ip)*dt 
p(ip)= gammal * ( rhoE(ip) - 0.5*rhov(ip)*v(ip) ) 
end if 
end do 
c 

c Output . 

c ======= 

c 

if ( mod(i, 10) .eq. 0 ) then 
write (50,*)’ it,dt = ’,i,dt 

write(50,*)’ it npoin(active) nelem(active) ’ 

c 

c How many active points (neac) ? 

c 

npac=0 

do ip=l,mpoin 

if ( ipact(ip) .eq. 1 ) npac=npac+l 
end do 
c 

c How many active elements (neac) ? 

c 

neac=0 

do ie=l,melem 

if ( mrhist(ie,5) .eq. 1 ) then 
neac=neac+l 

1 (intmat (ie , 1) )=mrhist (ie ,4) 
end if 
end do 

write (50,*) i,npac ,neac, npoin, nelem 
write(50,*)’ xp rho v 

& ’ rhoE p 1 ’ 

do ip=l, npoin 
if ( ipact(ip) .eq. 1 ) 

fe write(50,6) xp(ip) ,rho(ip) ,v(ip) ,rhoE(ip) ,p(ip) ,l(ip) 
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end do 
end if 
c 

1000 continue 


c 

c —END OF MAIN TIME STEP LOOP 

c 

5 f ormat (8(2x , f6 . 3) ) 

6 f ormat (5(2x ,el2 . 5) , lx, il) 
stop 

end 

c 


c 


c 


c 

c 

c 

c 

c 


subroutine init (gamma, cour.ntime ,di ,ctore ,ctode , epsil .nrmax , 
& nbuff ,ntref ,npoin, nelem, mpoin, melem, 

& period, intmat, mrhist, lrefe, ldere, ipact , 

& xp, rho, p, v, rhov, rhoE, 

& ibdry_l, ibdry_r ) 

real xp (mpoin) , rho (mpoin) ,p (mpoin) ,v(mpoin) , rhov (mpoin) 
real rhoE(mpoin) 

integer intmat (melem, 2) .mrhist (melem, 5) , lrefe (melem) 
integer ldere (melem) , ipact (mpoin) 


parameter (pi=3 . 141592654) 

Set constants, etc. 

gamma=l . 4 
gamma l=gamma-l . 
cour =0.5 
ntime=3000 
di =0.5 
ctore= . 25 
ctode= . 1 


Courant number 
Number of time steps 
Diffusion coefficient 
Threshold of error fen. to refine 
Threshold to de-refine 
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epsil=0 . 005 
nbuf f =3 
ntref =2 
nrmax=3 
c 

c Set the initial conditions 

c 

npoin=128 

nelem=127 

c 

do ie=l,npoin-l 
intmat (ie, l)=ie 
intmat (ie , 2) =ie+l 
end do 
c 

c boundary conditions. 

c 

if (period .ne. 0.) then 
intmat (nelem , 1) =nelem 
intmat (nelem , 2) =1 
else 

ibdry_l=l 
ibdry_r=npoin 
end if 
c 

do ie=l,melem 
do ir=l,5 
mrhist (ie , ir)=0 
end do 
lrefe(ie)=0 
ldere (ie)=0 
end do 
c 

do ie=l, nelem 
mrhist(ie,5)=l 
end do 
c 


Refinement filter coefficient 
No. of buffer layers 
No. of time steps per refinement 
Max. No. of Refinement Levels. 


16 



do ip=l,npoin 
ipact (ip) =1 
end do 
c 

do ip=npoin+l ,mpoin 
ipact(ip)=0 
end do 
c 

c Initial Configuration for Sod’s Shock Tube Problem 

c 

do ip=l,npoin/2 
rho(ip)= 8. 

p(ip)=10 . 

end do 
c 

do ip=npoin/2, npoin 
rho(ip)= 1. 
p(ip)= 1. 
end do 
c 

do ip=l,npoin 

xp(ip) =f loat (ip-1) *6 .4/float (npoin) 
v(ip) =0.0 
rhov(ip)=rho(ip)*v(ip) 

rhoE(ip)=p(ip)/gammal + 0.5*rho(ip)*v(ip)*v(ip) 
end do 
c 

c Write out the initial conditions 

c 

it=0 

write (50 , *) ’ initial conditions ’ 

write(50,*)’ it npoin nelem’ 

write(50,*) it ,npoin, nelem 

write(50,*) ’ xp rho v 

& ’ rhoE p ’ 

do i=l, npoin 

write(50,6) xp(i) ,rho(i) ,v(i) ,rhoE(i) ,p(i) 
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end do 


c 

6 

c 

c 

c 


f ormat (5(2x,el2 . 5) ) 
return 
end 


Q************************************************************** ************* 


c 

subroutine geom(mpoin, melem, npoin, nelem, intmat , xp, 


& ce, rmatm, rlen, period, mrhist, ipact ) 

c 

c COMPUTE THE GEOMETRY ARRAYS FOR A ID FE FLOW SIMULATION. 

c 


integer mpoin, melem, npoin, nelem, intmat (melem, 2) 
real xp(mpoin) ,ce(melem,2) , rmatm (mpoin) ,rlen(melem) 
integer mrhist (melem , 5) , ipact (mpoin) 
c 

c Initialize the mass matrix 

c 

do ip=l, npoin 
rmatm ( ip) =0. 
end do 


c 

c Initialize element arrays 

c 

do ie=l, nelem 
rlen(ie)=0 . 
end do 
c 

c Compute the length of the elements (rlen) 

c 

do ie=l, nelem 

if ( mrhist (ie, 5) .eq. 1 ) 

& rlen(ie)= xp ( intmat (ie, 2)) - xp(intmat(ie, 1) ) 
end do 
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c 


Periodic b.c.? 


c 
c 

if ( period .ne. 0. ) then 
do ie=l,nelem 

if ( abs(rlen(ie)) .gt. period/2. ) 

& rlen(ie)=period + rlen(ie) 

end do 
end if 
c 

c The basis functions’ gradient 

c 

do ie=l,nelem 

if ( mrhist(ie,5) .eq. 1 ) then 
ce(ie,2)=l ./rlen(ie) 
ce(ie, 1)=-1 . *ce(ie,2) 
end if 
end do 
c 

c The inverse lumped mass matrix. 

c 

do ie=l,nelem 

if ( mrhist(ie,5) .eq. 1 ) then 
do in=l,2 

rmatm(intmat (ie , in) ) = rmatm(intmat (ie , in) ) + rlen(ie) 
end do 
end if 
end do 
c 

c Invert the sum 

c 

do ip=l,npoin 
if ( ipact (ip) . eq. 1 ) 

& rmatm(ip) = 2./rmatm(ip) 
end do 
c 

return 
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end 


c 

c *********************************************************************** 

c 

subroutine grerror (nelem, npoin, melem, mpoin, vref , 


& tempea, tempec, temppa, temppb, temppc, 

& epsil, ctore, ctode, nbuff, errorp, 

& errore, mrhist, ipact, 

& lrefe, ldere, ce, rlen, intmat, 

& ietemp,iptemp,ibdry_l,ibdry_r ) 

c 

c THIS SUBROUTINE COMPUTES THE H-2 SEMINORM ERROR ESTIMATOR 

c [ errorp = temppa/ ( temppb + epsil*temppc) ] 

c (Ref : Lohner) 

c 


integer intmat (melem, 2) .mrhist (melem, 5) , ipact (mpoin) 
integer lrefe (melem) , ldere (melem) 
integer iptemp(mpoin) , ietemp (melem) 
real ce (melem , 2) , rlen (melem) 

real tempea(melem) , tempec (melem) .errore (melem) 

real temppa(mpoin) .temppb (mpoin) .temppc (mpoin) .errorp (mpoin) 

real vref (mpoin) 


c 

c nbuff =3 

c 

c Initialize element arrays 


c 

do ie=l, melem 
tempea(ie)=0 . 
tempec (ie)=0. 
errore(ie)=0 . 
lrefe(ie) =0. 
ldere (ie) =0. 
end do 
c 

c Initialize point arrays 

c 
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do ip=l,rapoin 
temppa(ip)=0 . 
temppb(ip)=0 . 
temppc(ip)=0. 
errorp(ip)=0 . 
end do 

-Sums at the elements . 
do ie=l,nelem 

if (mrhist(ie,5) .eq. 1) then 
do in=l,2 

tempea(ie) = tempea(ie) + ce(ie , in)*vref (intmat (ie , in) ) 
tempec(ie) = tempec(ie) + abs( ce(ie,in) ) * 

abs( vref (intmat (ie, in)) ) 

end do 
end if 
end do 

Scatter to the points, 
do ie=l,nelem 

if (mrhist(ie,5) .eq. 1) then 
do in=l,2 

temppa(intmat (ie , in) ) = temppa(intmat(ie,in)) + 

rlen(ie)*ce(ie,in)*tempea(ie) 
temppb ( intmat (ie , in) ) = temppb(intmat(ie,in)) + 

rlen(ie) * abs(ce(ie , in) ) * 
abs(tempea(ie)) 

temppc(intmat(ie,in)) = temppc(intmat(ie,in)) + 

rlen(ie) * abs(ce(ie , in) ) * 
tempec(ie) 

end do 
end if 
end do 


Boundary terms 



c 

temppa(ibdry_l) = 0. 
temppa(ibdry_r) = 0. 
c 

c The error estimate at the points 

c 

do ip=l,npoin 


if ( ipact(ip) .eq. 1 ) 

& errorp(ip) = abs (temppa(ip) ) / 

& ( temppb(ip) + epsil*temppc (ip) ) 

end do 
c 

c Compute the error estimate at the elements (max. of node values). 


c 

do ie=l,nelem 


if (mrhist(ie,5) .eq. 1) 

& errore(ie)=amaxl( errorp(intmat(ie, 1)) , 

& errorp(intmat(ie,2)) ) 

end do 
c 

c Mark the elements for refinement or de-refinement. 

c 


do ie=l,nelem 

if (mrhist(ie,5) .eq. 1) then 
if (errore(ie) .gt . ctore) lrefe(ie)=l 
if (errore (ie) . It . ctode) ldere(ie)=l 
end if 
end do 
c 

c Buffer Layers 

c 

do ib=l,nbuff 
do ip=l,npoin 
iptemp(ip)=0 
end do 
c 

do ie=l,nelem 
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if (mrhist(ie,5) .eq. 1) then 
do in=l,2 

if ( lrefe(ie) .eq. 1) 

& iptemp(intmat(ie,in)) = 1 

end do 
end if 
end do 
c 

do ie=l,nelem 

if (mrhist (ie , 5) . eq. 1) then 
do in=l,2 

if ( iptemp(intmat (ie , in) ) .eq. 1) 

& lrefe(ie) = 1 

end do 
end if 
end do 
end do 
c 

c Space Continuity on levels of refinement (only necessary on ID) 

c 

do iloop=l,5 
do ie=l,nelem 
ietemp(ie) = 0 
if (mrhist(ie,5) .eq. 1) 

& ietemp(ie) = mrhist(ie,4) + lrefe(ie) 
end do 
c 

do ip=l,npoin 
iptemp(ip)=0 
end do 
c 

do ie=l,nelem 

if (mrhist (ie, 5) .eq. 1) then 
do in=l,2 

if ( ietemp(ie) .gt. iptemp(intmat(ie,in)) ) 

& iptemp(intmat (ie , in) ) = ietemp(ie) 

end do 


23 



end if 
end do 

do ie=l,nelem 

if (mrhist(ie,5) .eq. 1) then 
do in=l,2 

if ( iptemp(intmat(ie,in))-mrhist(ie,4) .gt. 1) 
& lrefe(ie) = 1 

if ( iptemp(intmat(ie,in))-mrhist(ie,4) .ge. 1) 
& ldere(ie) = 0 

end do 

ietemp(ie) = mrhist(ie,4) + lrefe(ie) 
do in=l,2 

if ( ietemp(ie) .gt. iptemp(intmat(ie,in) ) ) 

& iptemp(intmat(ie,in)) = ietemp(ie) 

end do 

end if 
end do 
end do 

return 

end 


c 


c ***************** ************* ****** ****** ********************************* 
c 

subroutine unref(melem, mpoin, nelem, npoin, intmat. 


& mrhist, ldere.ipact ) 

c 

c THIS ROUTINE DE-REFINES THOSE ELEMENTS MARKED FOR DE-REFINEMENT 

c 


integer intmat (melem, 2) 
integer mrhist (melem, 5) , ldere (melem) 
integer ipact (mpoin) 
c 
c 
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c AMR arrays 

c ========== 

c 

c mrhist (ie , 1) == parent element number. (default=0) 

c 2) == first child element number. (default=0) 

c 3) == second child element number. (default=0) 

c 4) == level of refinement. (default=0) 

c 5) == element active? ( l=yes ) 

c 

c ldere(ie) =1 ======> de-refine 

c =0 ======> do not de-refine 

c 

c Loop over the elements 


c 

do ie=l,nelem 

if ( ldere (ie) . eq. 1 .and. mrhist(ie,4) .ne.O ) then 
c 

c Identify elements involved 

c 

iparent=mrhist (ie , 1) 
ichildl=mrhist (iparent , 2) 
ichild2=mrhist (iparent ,3) 


c 

c Derefine iff both children are marked for de-refinement (4 steps) 

c 

if ( ldere(ichildl) .eq. 1 .and. ldere (ichild2) . eq. 1 ) 

& then 

c 

c (1) Identify and ‘‘deactivate" points not assoc, w/parent. 


c 

do in=l,2 

ipl=intmat (ichildl , in) 
ip2=intmat ( ichild2 , in) 
if (ipl .ne . intmat (iparent , 1) .and. 

& ipl .ne . intmat (iparent ,2) ) ipact(ipl)=0 

if (ip2 .ne . intmat (iparent , 1) .and. 

& ip2 . ne . intmat (iparent ,2) ) ipact(ip2)=0 
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c 

c 


end do 


c 


c 

c 

c 

c 

c 

c 


..(2) Disable children elements... 

mrhist (ichildl ,5)=0 
mrhist (ichild2 ,5)=0 

..(3) Enable parent element 

mrhist (iparent , 5) =1 

. . (4) Unmark the children elements 

ldere (ichildl) =0 
ldere (ichild2) =0 


c 

c 

c 

c 


end if 

end if 

end do 

return 

end 


! ( ... if all children marked) 
! ( ... if one child marked) 

! (end loop over elements) 


c 

c************************************************************************ 

c 

subroutine ref ine(melem, mpoin, nelem, npoin, intmat, nrmax, 


& mrhist, xp, rho, rhov, rhoE, lrefe, ipact ) 

c 

C THIS ROUTINE REFINES THOSE ELEMENTS MARKED FOR REFINEMENT 


c 

integer intmat (melem, 2) 
real xp (mpoin) 

real rho (mpoin) , rhov (mpoin) .rhoE(mpoin) 
c 

integer mrhist (melem, 5) ,lrefe(melem) 
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integer ipact(mpoin) 
c 

c AMR arrays 

c ========== 

c 

c mrhist (ie , 1) == parent element number. (default=0) 

c 2) == first child element number. (default=0) 

c 3) == second child element number. (default=0) 

c 4) == level of refinement. (default=0) 

c 5) == element active? ( l=yes ) 

c 

c lrefe(ie) =1 ======> refine 

c =0 ======> do not refine 

c 

c Loop over the elements 


c 

nrmaxl=nrmax-l 
do ie=l,nelem 

if ( lrefe(ie) .eq. 1 .and. mrhist(ie,4) .le.nrmaxl ) then 
npoin=npoin+l !new point 

ipact(npoin)=l 
if (npoin .gt. mpoin) then 

write(*,*) ’Please increase mpoin. Needed: ’ ,npoin 
stop 
end if 
c 

nelem2=nelem+2 

if (nelem2 .gt. melem) then 

write (*,*) ’Please increase melem. Needed :’ ,nelem2 
stop 
end if 
c 

xp (npoin) = 0.5*( xp(intmat(ie,l)) + xp(intmat (ie , 2) ) ) 
rho(npoin) = 0.5*( rho(intmat(ie,l)) + 

& rho(intmat(ie,2) ) ) 

rhov(npoin) = 0.5*( rhov(intmat(ie,l)) + 

& rhov(intmat(ie,2)) ) 
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rhoE(npoin) = 0.5*( rhoE(intmat (ie , 1) ) + 

& rhoE(intmat(ie,2)) ) 

c 

c New Element Connectivity (Input new elements at the end of the list) 


c 

intmat (nelem+1 , 1) = intmat (ie,l) 
intmat (nelem+1, 2) = npoin 
intmat (nelem+2,1) = npoin 
intmat (nelem+2 ,2) = intmat (ie, 2) 
c 

c Disable the parent elements and enable the children 

c 

mrhist(ie,5)=0 
mrhist (nelem+1, 5) = 1 
mrhist (nelem+2, 5) = 1 
c 

c Modify the rest of the mesh refinement history parameters 

c 

mrhist (ie, 2) = nelem+1 
mrhist (ie, 3) = nelem+2 
mrhist (nelem+1 , 1) = ie 
mrhist (nelem+2 , 1) = ie 
mrhist (nelem+1 ,4) = mrhist (ie, 4) + 1 
mrhist (nelem+2, 4) = mrhist (ie, 4) + 1 
c 

c Modify the current number of elements 

c 

nelem=nelem+2 
end if 
end do 
c 

return 

end 

c 

c 

subroutine f luxld(mpoin, melem, npoin, nelem, gamma, cour, 


28 



& intmat , rho, rhov, rhoE, p, v, 

& rlen, vraax, flux, dt, mrhist, ipact ) 

c 

c Compute the fluxes at the points and the time step. 

c =================================================== 

c 


integer mpoin , melem , npoin , nelem , intmat (melem , 2) 
real rho(mpoin) ,rhov(mpoin) ,rhoE(mpoin) 
real p(mpoin) ,v(mpoin) ,f lux(mpoin,3) 
real rlen (melem) ,vmax (mpoin) 
integer mrhist (melem, 5) , ipact (mpoin) 
c 

cbig=l . e30 
gamma l=gamma-l . 
do ip=l, npoin 
if ( ipact (ip) .eq. 1 ) then 
v(ip)=rhov(ip)/rho(ip) 

p(ip)= gammal * ( rhoE(ip) - 0 . 5*rhov(ip) *v(ip) ) 
c 

flux(ip,l) = rhov(ip) 
flux(ip,2) = rhov(ip)*v(ip) + p(ip) 
flux(ip,3) = rhoE(ip) *v(ip) + p(ip) *v(ip) 
end if 
end do 


c 

c TIME STEP: dt . 

c ============== 

c 

c Compute the ‘‘maximum velocity" (sound speed + v) at the points. 

c 


do ip=l, npoin 
if ( ipact (ip) . eq. 1 ) 

& vmax(ip) = sqrt (gamma*p(ip) /rho (ip) ) + abs(v(ip)) 
end do 
c 

c and the time step. . . 

c 


29 



dt= cbig 
do ie=l,nelem 

if ( mrhist (ie, 5) .eq. 1 ) then 
do in=l,2 

dt = aminl( dt, rlen(ie)/vmax(intmat(ie,in)) ) 
end do 
end if 
end do 
c 

dt = cour * dt 
c 

return 

end 

c 

c**************************** *********** *********************************** 


c 

subroutine deltau(mpoin, melem, npoin.nelem, intmat , 

& du , dup .flux , dt , ce , rlen , cour , di , rho , rhov , rhoE , dif f , 

ft vmax, mrhist) 

c 

c Compute the du ’ s and dup ’ s 


c 

integer mpoin, melem, npoin.nelem, intmat (melem, 2) 
real du(melem,3), dup (mpoin, 3) .flux (mpoin, 3) ,dt 
real ce (melem, 2) , rlen (melem) 

real rho (mpoin) , rhov (mpoin) .rhoE(mpoin) , dif f (melem, 3) 
real vmax (mpoin) 
integer mrhist (melem, 5) 
c 

cbig=l .e20 
csmall=-l .e20 
c 

c Initialize du 

c 

do im=l,3 
do ie=l,nelem 
du(ie,im) = 0 
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diff(ie,im) = 0. 
end do 
end do 
c 

c The minimum delta x. 

c 

rlmin=cbig 
rlmax=c small 
do ie=l,nelem 

if (mrhist(ie,5) .eq. 1 .and. rlen(ie) . It .rlmin) 

& rlmin=rlen(ie) 

if (mrhist(ie,5) .eq. 1 .and. rlen(ie) .gt .rlmax) 

& rlmax=rlen(ie) 

end do 
c 

c Initialize dup 

c 

do im=l,3 
do ip=l,npoin 
dup (ip , im)=0 . 
end do 
end do 
c 

c BEGUN LOOP OVER THE ELEMENTS 

c 

do 1000 ie=l,nelem 
if ( mrhist(ie,5) .eq. 1 ) then 
c 

c Compute du at the elements 

c 

do im=l,3 
do in=l,2 

du(ie.im) = du(ie,im) + ce(ie, in)*f lux(intmat (ie , in) , lm) 
end do 
end do 
c 

c Compute the diffusion terms at the elements 
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c 

vv=amaxl (vmax(intmat (ie , 1) ) , vmaxdntmat (ie , 2) ) ) 
do in=l,2 

diff(ie,l) = diff(ie,l) + ce(ie,in)*rho(intmat(ie,in)) 

& *vv 

diff(ie,2) = diff(ie,2) + ce(ie,in)*rhov(intmat (ie , in) ) 
& *vv 

diff(ie,3) = diff(ie,3) + cede ,in)*rhoE(intmat (ie , in) ) 
& *vv 

end do 
c 

dfact = 0.5*di 
c 

do im=l,3 

diff(ie,ini) = dfact * dif f (ie , im) *rlen(ie) *rlmin 
end do 
c 

c Scatter to the points. 

c 

do im=l,3 
do in=l,2 

dup(intmat(ie,in) ,im) = dupdntmat (ie,in) , im) + 

& 0.5*rlen(ie)*du(ie,im) + 

& ce(ie,in)*diff (ie,im) 

end do 
end do 
c 

end if 
c 

1000 continue 
c 

c END OF THE ELEMENT LOOP 

c 

return 

end 
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